Load required libraries
library(tidyverse)
── Attaching core tidyverse packages ───────────────────────────────────────────────────────────────────────────────────────── tidyverse 2.0.0 ──
✔ dplyr 1.1.4 ✔ readr 2.1.5
✔ forcats 1.0.0 ✔ stringr 1.5.1
✔ ggplot2 3.5.2 ✔ tibble 3.2.1
✔ lubridate 1.9.4 ✔ tidyr 1.3.1
✔ purrr 1.0.4
── Conflicts ─────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag() masks stats::lag()
ℹ Use the ]8;;http://conflicted.r-lib.org/conflicted package]8;; to force all conflicts to become errors
library(here)
G3;here() starts at /Users/peter.macpherson/Library/Mobile Documents/com~apple~CloudDocs/Documents/Documents - Peter’s Mac mini/Projects/AQ-MTB_2025-05-14/Work/Data/aq_mtb
g
library(brms)
G3;Loading required package: Rcpp
gG3;Loading 'brms' package (version 2.22.0). Useful instructions
can be found by typing help('brms'). A more detailed introduction
to the package is available through vignette('brms_overview').
gG3;
Attaching package: ‘brms’
gG3;The following object is masked from ‘package:stats’:
ar
g
library(tidybayes)
G3;
Attaching package: ‘tidybayes’
gG3;The following objects are masked from ‘package:brms’:
dstudent_t, pstudent_t, qstudent_t, rstudent_t
g
library(sf)
G3;Linking to GEOS 3.13.0, GDAL 3.8.5, PROJ 9.5.1; sf_use_s2() is TRUE
g
library(mgcv)
G3;Loading required package: nlme
gG3;
Attaching package: ‘nlme’
gG3;The following object is masked from ‘package:dplyr’:
collapse
gG3;This is mgcv 1.9-3. For overview type 'help("mgcv-package")'.
gG3;
Attaching package: ‘mgcv’
gG3;The following objects are masked from ‘package:brms’:
s, t2
g
library(priorsense)
library(osmdata)
G3;Data (c) OpenStreetMap contributors, ODbL 1.0. https://www.openstreetmap.org/copyright
g
library(nngeo)
G3;Registered S3 method overwritten by 'data.table':
method from
print.data.table
g
Import the modelling dataset - cleaned and prepped in the
2025-05-12_aq_clean.Rmd script
aq_model_data <- read_rds("input_data/aq_model_data.rds")
#also load the scale clusters
load("input_data/scale_72_clusters.rda") #SCALE clusters
#and the grid data
read_rds("input_data/mw_100m_cropped.rds")
$mwi_ppp_2020.tif
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,] 0.045165006 0.046947237 0.055178735 0.047102101 0.043815471 0.041929662 0.039621390 0.03764743 3.824615e-02 3.529998e-02
[2,] 0.054517046 0.047781128 0.052100033 0.046856463 0.042057369 0.028544975 0.027005022 0.02849834 2.889967e-02 2.749119e-02
[3,] 0.045086842 0.047777772 0.041951492 0.034255959 0.029440563 0.024227845 0.016665278 0.02430716 7.012391e-03 1.374129e-02
[4,] 0.035714723 0.039766397 0.035008151 0.032815386 0.030017985 0.028784296 0.029581040 0.02444118 1.111266e-02 4.598498e-03
[5,] 0.015322668 0.018739052 0.018032806 0.018089430 0.017559310 0.018649841 0.017364351 0.01674390 1.384361e-02 5.014876e-03
[,11] [,12] [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20]
[1,] 3.581402e-02 0.033903446 0.028697260 0.025856392 9.365267e-03 4.969447e-03 4.880766e-03 6.018387e-03 7.693967e-03 1.687009e-02
[2,] 3.080107e-02 0.031454172 0.036637034 0.028348375 1.773694e-02 5.785842e-03 6.604420e-03 6.460932e-03 8.471724e-03 1.088048e-02
[3,] 2.850383e-02 0.033486295 0.035977926 0.021789778 7.277420e-03 5.582869e-03 6.420635e-03 6.296951e-03 8.607163e-03 7.565869e-03
[4,] 4.251544e-03 0.014496404 0.024309691 0.025809046 5.967221e-03 5.781442e-03 5.793448e-03 8.620093e-03 8.066839e-03 7.239441e-03
[5,] 3.449674e-03 0.004229320 0.027064180 0.021735243 2.350955e-02 6.558061e-03 5.777395e-03 6.241026e-03 6.318555e-03 7.906212e-03
[,21] [,22] [,23] [,24] [,25] [,26] [,27] [,28] [,29] [,30]
[1,] 7.483703e-03 1.134669e-02 2.477654e-02 1.747209e-02 0.02301022 0.03323989 0.03443765 0.04373282 0.03239537 0.03268054
[2,] 6.492021e-03 1.798505e-02 3.820698e-02 3.622469e-02 0.04873810 0.06264146 0.04502079 0.04902079 0.03930074 0.04743655
[3,] 1.413537e-02 2.868031e-02 4.856283e-02 7.635828e-02 0.10526711 0.10480442 0.09942765 0.08063150 0.07026784 0.07916132
[4,] 3.765117e-02 4.892825e-02 8.083706e-02 6.265849e-02 0.06352694 0.10991651 0.11036292 0.10957718 0.09102185 0.09063222
[5,] 4.181869e-02 6.396030e-02 6.208089e-02 5.934627e-02 0.06683331 0.11021938 0.11196631 0.10069560 0.08459925 0.09193873
[,31] [,32] [,33] [,34] [,35] [,36] [,37] [,38] [,39] [,40]
[1,] 0.01216105 6.328823e-03 0.01512775 0.02006114 3.101150e-02 0.05473484 0.06160709 0.07647948 0.07574075 0.07371373
[2,] 0.04709647 1.577023e-02 0.01547693 0.01512168 9.752575e-03 0.03868737 0.05658274 0.06831920 0.07236640 0.06479138
[3,] 0.08027354 4.902577e-02 0.04347426 0.04361450 2.506110e-02 0.02168380 0.03363772 0.07003722 0.10318663 0.07466165
[4,] 0.10364193 1.078961e-01 0.06888389 0.06897457 7.109120e-02 0.06687172 0.04877206 0.06761945 0.06558794 0.05342549
[5,] 0.10197861 1.167887e-01 0.11340023 0.07957450 8.267088e-02 0.07813363 0.05316366 0.03048283 0.04607106 0.05050959
[,41] [,42] [,43] [,44] [,45] [,46] [,47] [,48] [,49] [,50]
[1,] 0.09937067 0.09119225 0.08600094 0.09103850 0.09306126 0.17105776 0.16657399 0.19017223 0.18780756 0.25865054
[2,] 0.08613592 0.09276834 0.08768278 0.09463739 0.09931904 0.20142318 0.34862801 0.24239205 0.20016080 0.14733949
[3,] 0.08796223 0.08529528 0.05861057 0.08760586 0.17597923 0.18829565 0.30547777 0.18160681 0.10488693 0.05029128
[4,] 0.07356653 0.04595139 0.06095412 0.07659898 0.20957065 0.22002560 0.26929030 0.13957511 0.13037296 0.12920222
[5,] 0.05924074 0.12837498 0.10924210 0.04965207 0.03692738 0.13063972 0.13985389 0.22550066 0.22139461 0.23821150
[,51] [,52] [,53] [,54] [,55] [,56] [,57] [,58] [,59] [,60]
[1,] 0.18521543 1.718710e-01 0.17261122 1.592420e-01 0.13821790 2.070713e-01 2.180091e-01 2.510735e-01 3.347855e-01 0.42194489
[2,] 0.25576282 1.704141e-01 0.16087675 6.019334e-02 0.05911760 2.086819e-01 2.241761e-01 2.393655e-01 2.442953e-01 0.28544879
[3,] 0.08434317 1.395781e-01 0.03054703 2.303165e-02 0.02988590 1.481759e-01 1.958029e-01 2.119675e-01 2.289245e-01 0.22243160
[4,] 0.12479196 1.226818e-01 0.12564214 1.318173e-01 0.13020913 1.327822e-01 2.432093e-01 2.560078e-01 2.778168e-01 0.30719724
[5,] 0.22799887 1.374678e-01 0.14205216 1.374896e-01 0.14270458 2.208541e-01 2.599400e-01 2.764551e-01 4.585261e-01 0.61238015
[,61] [,62] [,63] [,64] [,65] [,66] [,67] [,68] [,69] [,70] [,71]
[1,] 0.20092252 0.2094545 0.23813984 0.2600536 0.28496569 0.3515317 0.3495567 3.0187855 3.4292254 4.3702903 10.2430563
[2,] 0.37174022 0.2000832 0.20835581 0.2420900 0.25386712 0.3162111 0.3505953 2.8476219 2.6364815 3.5290518 4.9171147
[3,] 0.31947485 0.3825941 0.37987089 0.2292890 0.24251479 0.2557049 2.4554725 2.3886611 1.5429624 2.8309567 3.7996442
[4,] 0.34511879 0.5266556 0.32725015 0.2879633 0.32016167 0.8491008 2.2822711 3.1561983 4.8729448 4.3276281 5.2610259
[5,] 1.26680827 1.3391285 1.04976392 0.3362455 0.29695830 0.2514963 0.7431583 6.0843654 6.2483072 5.5046158 5.5970216
[,72] [,73] [,74] [,75] [,76] [,77] [,78] [,79] [,80] [,81] [,82]
[1,] 12.3174086 20.4170475 5.9770074 5.8101201 6.1680431 7.0288363 2.9920485 2.7822750 3.2919111 2.7028601 2.7344890
[2,] 11.4855623 3.7609022 3.5631394 2.8713186 2.5883286 6.6772809 6.6999898 6.8617449 3.1268249 2.7026587 2.1398206
[3,] 4.0249171 2.3710566 2.3969679 1.5280255 1.3292254 1.3552140 6.3890018 6.8920774 5.9485264 6.3546324 2.1996837
[4,] 1.3237036 1.1526742 2.2051387 1.6958532 5.7821097 1.2228683 2.1854522 6.9267850 5.9227381 5.7966475 2.1904876
[5,] 1.9726193 1.2746783 1.6196227 1.4951222 5.7002840 1.2784475 5.3949475 2.2073100 2.0879290 2.0115442 1.1903656
[,83] [,84] [,85] [,86] [,87] [,88] [,89] [,90] [,91] [,92] [,93]
[1,] 2.8262646 3.1302779 3.6975005 2.1694636 2.2199211 2.3952386 2.0133886 1.0360730 6.0406532 6.6095033 6.3419876
[2,] 1.7706543 1.8591148 2.6041703 1.9700989 2.0043545 2.1224627 12.0202541 10.6813793 5.7113476 5.8507247 5.5529652
[3,] 0.7009542 0.5428279 1.0550138 0.9149451 9.6139355 10.1034241 9.7381792 7.6451182 5.3572125 5.1382504 5.1880946
[4,] 0.6239597 0.5024059 0.9726501 0.9038598 9.6681395 7.2640862 6.9079909 5.7228518 6.3007116 6.8418317 7.1582618
[5,] 0.4369978 0.4899155 0.8078130 0.7381263 8.1139145 6.2780352 5.4850063 5.4929547 6.6182308 8.1480541 7.9830341
[,94] [,95] [,96] [,97] [,98] [,99] [,100] [,101] [,102] [,103] [,104] [,105]
[1,] 6.0880580 4.9141769 5.6921430 5.6045480 5.6583200 5.7986665 4.8305206 4.5144963 4.3507361 4.8232894 4.1951013 4.6705604
[2,] 5.3891821 5.2441683 4.3157148 4.7750669 4.9417787 5.7972479 5.7227149 6.0235286 5.1872697 5.3412123 5.1507735 4.7766061
[3,] 5.2370167 5.1144838 4.4092488 4.2801752 4.8011594 5.8726711 5.7670879 6.0404143 6.7926741 6.6716657 5.2236753 5.0849514
[4,] 7.2552085 6.9078765 5.8223286 5.2468567 5.3122015 6.6020727 6.6301870 6.3754029 7.2572618 7.1930847 7.2747545 7.0348992
[5,] 7.6431541 7.6273088 6.5719180 5.9137011 5.9537210 6.5491629 6.6482563 6.1040545 6.8009992 7.0888548 7.2325339 6.9409246
[,106] [,107] [,108] [,109] [,110] [,111] [,112] [,113] [,114] [,115] [,116]
[1,] 4.0513868 4.5824718 31.0151443 35.9615288 31.9965153 33.5947914 34.8950462 33.6029587 32.1246338 31.0912933 13.1747513
[2,] 4.3982000 4.5165482 30.4003181 31.6037254 30.7493134 31.7781143 35.9739647 37.8679657 37.2548790 15.6511335 14.6784821
[3,] 5.1197200 4.4855256 33.6331902 39.4704590 25.6988716 33.0550232 13.5869970 13.1186905 13.0235996 13.0852079 16.8719463
[4,] 7.5529518 6.8021359 38.6025314 38.8714218 16.4733143 13.1496153 11.0804739 9.1725912 11.2854414 11.4492359 14.8278837
[5,] 7.0277257 7.7998891 6.9491796 14.6560345 17.5421658 16.1886940 11.8365650 11.4074469 13.3405104 13.5668907 14.1344957
[,117] [,118] [,119] [,120] [,121] [,122] [,123] [,124] [,125] [,126] [,127]
[1,] 17.4172173 14.5559740 10.0551395 0.9203696 0.4967222 0.4107275 0.3978702 0.4675624 0.7120231 0.5725706 0.5796408
[2,] 16.7342072 11.8646288 11.6098852 0.8615871 0.6930147 0.6287227 0.4980069 0.7469794 0.5880303 0.4945861 0.4778275
[3,] 17.6651497 16.5033340 14.2913799 1.1762685 1.0414619 0.9750723 0.9179240 0.6648465 0.6006990 0.4671351 0.4315493
[4,] 15.0737171 14.5026913 1.4341162 1.7666528 1.2385327 1.2163831 1.1903278 1.0431690 0.6815187 0.6674115 0.6116236
[5,] 14.0739212 0.9036093 1.0228646 1.7256023 1.1906258 1.1074116 1.0897969 0.8469090 0.8302435 0.7722517 0.7019941
[,128] [,129] [,130] [,131] [,132] [,133] [,134] [,135] [,136] [,137] [,138] [,139]
[1,] 0.5407016 5.6824121 5.5070438 5.393878 5.314744 5.378588 5.516953 5.472218 6.871101 7.120472 2.490248 1.913114
[2,] 0.4782741 0.5772134 5.7760925 5.619320 5.583245 5.598059 5.671175 6.422924 7.114538 6.954295 2.615454 2.365523
[3,] 0.4778309 6.0533495 5.8488727 5.827873 5.702469 5.692585 2.055359 2.058331 2.605898 2.534510 2.526598 1.925066
[4,] 6.0364213 6.8703146 6.9667211 6.728249 2.507176 2.410894 2.355129 2.363863 2.335734 2.941173 2.280494 2.165629
[5,] 7.4033751 6.2904220 6.0554729 5.708776 2.480082 2.428432 2.440901 2.502323 2.452162 2.456221 2.371370 2.334642
[,140] [,141] [,142] [,143] [,144] [,145] [,146] [,147] [,148] [,149] [,150] [,151]
[1,] 1.890159 1.864571 1.824639 1.813932 1.855493 1.862693 1.844257 1.884519 1.856309 1.911597 1.900523 1.8835605
[2,] 1.930770 1.902335 1.879635 1.835053 1.837441 1.848145 1.740682 1.844983 1.900461 1.933071 1.847667 1.8871386
[3,] 1.854827 1.917260 1.922128 1.888102 1.913544 1.864133 1.771954 1.868454 1.890933 1.831353 1.933270 1.9211165
[4,] 2.187655 2.135860 2.126067 2.083079 2.073362 2.087825 2.053761 2.042283 2.012455 2.036653 2.052104 1.9967043
[5,] 2.811734 2.288969 2.239703 2.176805 2.207034 2.119645 2.070684 2.082571 1.989093 2.088601 2.078266 2.0321856
[,152] [,153] [,154] [,155] [,156] [,157] [,158] [,159] [,160] [,161] [,162] [,163]
[1,] 1.8760191 4.2506528 1.8467858 1.9328438 4.556878 4.581502 5.842960 5.753704 4.339655 5.3159437 3.6693192 3.5465922
[2,] 1.8791692 4.2897072 4.4196887 4.4844222 4.750791 5.543705 5.940609 5.749850 4.410793 5.4709144 5.1160870 3.6569257
[3,] 1.8962737 4.3825073 4.1637969 4.6726365 4.731542 4.696297 5.871382 5.844340 4.788056 5.8306942 3.9437084 3.8018162
[4,] 2.6364172 2.5760887 4.5244684 4.6207194 5.156001 5.343033 5.413106 5.551797 5.356392 3.1116621 3.1301842 3.1519930
[5,] 2.6615610 2.5835884 2.7367139 2.7678025 5.207396 5.004921 4.730000 4.490225 6.635791 3.2056251 3.5913785 3.6644924
[,164] [,165] [,166] [,167] [,168] [,169] [,170] [,171] [,172] [,173] [,174] [,175] [,176]
[1,] 3.4214435 3.2192101 3.0704648 3.2135775 4.2359557 4.4121046 4.179620 5.511585 7.292460 4.900838 4.766211 4.537093 3.3284004
[2,] 3.6246648 3.5184062 3.2872734 3.1027408 4.3657980 5.0197420 6.747860 8.048831 7.825340 7.773330 4.906329 4.761338 3.9066331
[3,] 3.5798521 3.6366291 3.9564095 5.8185658 6.5702181 6.7610564 6.643492 7.909804 7.872230 8.015193 5.043985 5.217367 6.0537314
[4,] 2.3774395 3.6038818 3.4744658 4.0604777 4.2596602 4.2969036 4.415190 5.865735 7.650286 7.599905 7.653657 5.146161 5.9680948
[5,] 3.3508301 3.3644121 3.6195219 4.1825180 4.3101077 4.4454837 5.979053 9.620084 15.794232 10.375804 9.111597 5.311432 5.9213576
[,177]
[1,] 3.8432529
[2,] 9.6162081
[3,] 11.1199265
[4,] 11.5795832
[5,] 11.8466797
[ reached 'max' / getOption("max.print") -- omitted 172 rows ]
attr(,"dimensions")
$x
$from
[1] 2755
$to
[1] 2931
$offset
[1] 32.67125
$delta
[1] 0.0008333333
$refsys
Coordinate Reference System:
User input: WGS 84
wkt:
GEOGCRS["WGS 84",
ENSEMBLE["World Geodetic System 1984 ensemble",
MEMBER["World Geodetic System 1984 (Transit)"],
MEMBER["World Geodetic System 1984 (G730)"],
MEMBER["World Geodetic System 1984 (G873)"],
MEMBER["World Geodetic System 1984 (G1150)"],
MEMBER["World Geodetic System 1984 (G1674)"],
MEMBER["World Geodetic System 1984 (G1762)"],
MEMBER["World Geodetic System 1984 (G2139)"],
MEMBER["World Geodetic System 1984 (G2296)"],
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]],
ENSEMBLEACCURACY[2.0]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
CS[ellipsoidal,2],
AXIS["geodetic latitude (Lat)",north,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433]],
AXIS["geodetic longitude (Lon)",east,
ORDER[2],
ANGLEUNIT["degree",0.0174532925199433]],
USAGE[
SCOPE["Horizontal component of 3D system."],
AREA["World."],
BBOX[-90,-180,90,180]],
ID["EPSG",4326]]
$point
[1] FALSE
$values
NULL
attr(,"class")
[1] "dimension"
$y
$from
[1] 7627
$to
[1] 7803
$offset
[1] -9.363749
$delta
[1] -0.0008333333
$refsys
Coordinate Reference System:
User input: WGS 84
wkt:
GEOGCRS["WGS 84",
ENSEMBLE["World Geodetic System 1984 ensemble",
MEMBER["World Geodetic System 1984 (Transit)"],
MEMBER["World Geodetic System 1984 (G730)"],
MEMBER["World Geodetic System 1984 (G873)"],
MEMBER["World Geodetic System 1984 (G1150)"],
MEMBER["World Geodetic System 1984 (G1674)"],
MEMBER["World Geodetic System 1984 (G1762)"],
MEMBER["World Geodetic System 1984 (G2139)"],
MEMBER["World Geodetic System 1984 (G2296)"],
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]],
ENSEMBLEACCURACY[2.0]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
CS[ellipsoidal,2],
AXIS["geodetic latitude (Lat)",north,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433]],
AXIS["geodetic longitude (Lon)",east,
ORDER[2],
ANGLEUNIT["degree",0.0174532925199433]],
USAGE[
SCOPE["Horizontal component of 3D system."],
AREA["World."],
BBOX[-90,-180,90,180]],
ID["EPSG",4326]]
$point
[1] FALSE
$values
NULL
attr(,"class")
[1] "dimension"
attr(,"raster")
$affine
[1] 0 0
$dimensions
[1] "x" "y"
$curvilinear
[1] FALSE
$blocksizes
x y
[1,] 3892 1
attr(,"class")
[1] "stars_raster"
attr(,"class")
[1] "dimensions"
attr(,"class")
[1] "stars"
“We fitted a spatially smooth regression model for log-transformed PM2.5 concentrations using a Gaussian process (GP) smooth over spatial coordinates (x, y) to capture flexible spatial trends. Seasonal variation was modelled using a Fourier series expansion of day-of-year up to the 3rd harmonic (i.e., sine and cosine terms for annual, semi-annual, and tri-annual cycles) to account for complex seasonal patterns. The model was fitted in a Bayesian framework using brms with a Gaussian likelihood and the cmdstanr backend.”
Note, priors are still causing issues for me, so we will comment out now, and use the defauly stan weakly informative priors
Priors
priors <- c(
prior(normal(3.4, 1), class = "Intercept"),
prior(normal(0,1), class = "b"),
prior(exponential(1), class = "sigma"),
prior(exponential(1), class = "sdgp"),
prior(normal(0, 1), class = "lscale", lb = 0)
)
Fit prior only model
Now model with data.
pp_check(m1)
G3;Using 10 posterior draws for ppc type 'dens_overlay' by default.
g
Compare prior-only model to model with data using
priorsense package:
Predictions by space and time. Here, although we only have data within clusters, we will predict for cells outside of clusters based on covariates. We will also predict for the month for whcih we do not have data.
#set_up prediction date frame
nd_m1 <- st_as_sf(mw_100m_cropped) %>%
mutate(centroid = st_centroid(geometry)) %>%
mutate(x = st_coordinates(centroid)[, 1],
y = st_coordinates(centroid)[, 2]) %>%
select(-centroid) %>%
st_drop_geometry() %>%
rename(pop_density = 1) %>%
mutate(pop_density_km2 = pop_density / 0.01)
G1;H1;Errorh: object 'mw_100m_cropped' not found
Error during wrapup: not that many frames on the stack
Error: no more error handlers available (recursive errors?); invoking 'abort' restart
g
Now just draw 50 random coordinates, and predict by time to see whether we captured the trend.
#sample coordinate points from the household dataset
#we will fix the small number sampled later!
set.seed(123)
sampled_points_m1 <- aq_model_data %>%
sample_n(50) %>%
select(x, y)
#generate DOY 0–365 and calculate Fourier terms
doy_grid <- tibble(doy = 0:365) %>%
mutate(
sin_doy1 = sin(2 * pi * doy / 365),
cos_doy1 = cos(2 * pi * doy / 365),
sin_doy2 = sin(4 * pi * doy / 365),
cos_doy2 = cos(4 * pi * doy / 365),
sin_doy3 = sin(6 * pi * doy / 365),
cos_doy3 = cos(6 * pi * doy / 365)
)
#expand grid across sampled locations
prediction_df_m1 <- sampled_points_m1 %>%
crossing(doy_grid)
#add predictions
preds_m1 <- add_epred_draws(object = m1, newdata = prediction_df_m1)
#summarise
preds_m1_sum <- preds_m1 %>%
ungroup() %>%
group_by(doy) %>%
summarise(.epred_mean = mean(.epred),
.lower = quantile(.epred, probs=0.025),
.upper = quantile(.epred, probs = 0.975))
#GEt the observed data
observed_data <- aq_model_data %>%
select(mean_doy, log_pm2_5)
#plot
preds_m1_sum %>%
ggplot() +
geom_ribbon(aes(x=doy, ymin=.lower, ymax = .upper), fill="steelblue", alpha=0.3) +
geom_line(aes(x=doy, y=.epred_mean)) +
geom_jitter(data = observed_data, aes(x = mean_doy, y = log_pm2_5),
color = "darkred", alpha = 0.6, width = 0.5, height = 0.0, size = 1.2) +
labs(title = "Model-estimated log(PM2.5) with empirical measurements",
subtitle = "Model predictions with 95% CrI and observed data points",
x = "Day of year",
y = "log(PM2.5)",
caption = "Modelled estimates restricted to within clusters") +
theme_ggdist() +
theme(panel.background = element_rect(colour = "grey78"),
legend.position = "none")
Now we include grid level covariates of distance to the road, population density, and building footprint percent, along witht the mean temp and huidity on the day of measurement (from the purple air monitors)
Again, priors to be fixed later
Predictions by space and time
Sample coordinates, and plot over time
#sample coordinate points
#will sort out the small number later
set.seed(123)
sampled_points_m2 <- aq_model_data %>%
sample_n(50) %>%
select(x, y, mean_temp_c, mean_current_humidity, pop_density_km2, building_coverage_pct, dist_to_road_m, grid_id, log_pm2_5)
#generate DOY 0–365 and calculate Fourier terms
doy_grid <- tibble(doy = seq(0,365, by=1)) %>%
mutate(
sin_doy1 = sin(2 * pi * doy / 365),
cos_doy1 = cos(2 * pi * doy / 365),
sin_doy2 = sin(4 * pi * doy / 365),
cos_doy2 = cos(4 * pi * doy / 365),
sin_doy3 = sin(6 * pi * doy / 365),
cos_doy3 = cos(6 * pi * doy / 365)
)
#expand grid across sampled locations
prediction_df_m2 <- sampled_points_m2 %>%
crossing(doy_grid)
#add predictions
preds_m2 <- add_epred_draws(object = m2, newdata = prediction_df_m2, ndraws = 100)
#observed data for plotting
observed_data <- aq_model_data %>%
select(mean_doy, log_pm2_5) %>%
mutate(mean_doy = round(mean_doy))
#plot
preds_m2 %>%
ungroup() %>%
ggplot(aes(x=doy)) +
stat_lineribbon(aes(y=.epred), .width=0.95) +
geom_jitter(data = observed_data, aes(x = mean_doy, y = log_pm2_5),
color = "darkred", alpha = 0.6, width = 0.5, height = 0.0, size = 1.2) +
scale_fill_brewer() +
labs(title = "Model-estimated log(PM2.5) with empirical measurements",
subtitle = "Model predictions with 95% CrI and observed data points",
x = "Day of year",
y = "log(PM2.5)",
caption = "Modelled estimates restricted to within clusters") +
theme_ggdist() +
theme(panel.background = element_rect(colour = "grey78"),
legend.position = "none")
Compare models, using LOO CV
library(loo)
loo_m1 <- loo(m1)
loo_m2 <- loo(m2)
loo_compare(loo_m1, loo_m2)
TODO